

/*
program to calculate preferred habitat regressions
alternative specification

run_ph_alt_regression
inputs: 
- first variable (dep): change in yields D_shock
- second variable (ind): proxy for demand shock B_shock

- mat_var: maturity remaining (years)
- ra_var: risk aversion high/low indicator
-- note: ra_var may be (semi-)cts: 0 in non-crisis, but cts in crisis
- long_var: term length long/short indicator
- security_id: variable with security identifier (eg cusip)
- auction_id: variable with auction identifier (eg datetime)

- cluster_type: T or NT for clustering on auction or two-way auction/security
-- optional: include bw(#) for HAC-robust SEs

- mat_window: size of maturity window to include in each regression (years)
- mat_long_window (optional): size of maturity window for maturities>20 
- mat_step: size of step forward while rolling over maturities (years)

- fname_stub (optional): save results in rolling data directory

outputs:
- results_low, results_high: matrices with regression output by maturity


*/

********************************************************************************
********************************************************************************
cap program drop run_ph_alt_regression 
program define run_ph_alt_regression

* parameter arguments, maturity, and forward/zero rate variable names
syntax varlist , ///
	mat_var(varname) ra_var(varname) long_var(varname) ///
	security_id(varname) auction_id(varname) ///
	cluster_type(string) mat_window(real) mat_step(real) ///
	[mat_long_window(real 0) bw(real 0) fname_stub(string)]

***********************************
* first input variable: change in yields D_shock
local D_shock : word 1 of `varlist'
* second input variable: proxy for demand shock B_shock
local B_shock : word 2 of `varlist'

* shock interactions (allowing for cts ra_var)
sort `auction_id' `mat_var' `security_id'
gen B_shock_S0 = `B_shock'*(`long_var'==0)*(`ra_var'==0)
gen B_shock_S1 = `B_shock'*(`long_var'==0)*(`ra_var')
gen B_shock_L0 = `B_shock'*(`long_var'==1)*(`ra_var'==0)
gen B_shock_L1 = `B_shock'*(`long_var'==1)*(`ra_var')

* create panel id: group security by short/long and low/high
tempvar panel_id
egen `panel_id' = group(`ra_var' `long_var' `security_id')

* maximum maturity included in regressions
summ `mat_var'
local mat_max `r(max)'
* for "rolling" windows: optionally increase for >20y maturities
if `mat_long_window' == 0 {
	local mat_long_window = `mat_window'
}
* final left point of maturity windows
local mat_window_final = `mat_max' - `mat_long_window'

***********************************
* variance-covariance options
* note: if bw specified, must use ivreghdfe (slower)
* xtset data (note: only matters when bw option specified)
tempvar T
sort `auction_id' `mat_var' `security_id'
egen `T' = group(`auction_id')
xtset `panel_id' `T'

if `bw'==0 {
	local reg_cmd reghdfe
	if "`cluster_type'"=="T" {
		* clustered on auction
		local vce_opts vce(cluster `T')
	}
	else if "`cluster_type'"=="NT" {
		* two-way clustered on auction/security
		local vce_opts vce(cluster `T' `panel_id')
	}
	else if "`cluster_type'"=="r" {
		* basic robust SEs
		local vce_opts vce(robust)
	}
}
else {
	local reg_cmd ivreghdfe
	if "`cluster_type'"=="T" {
		* clustered on auction
		local vce_opts cluster(`T') bw(`bw')
	}
	else if "`cluster_type'"=="NT" {
		* two-way clustered on auction/security
		local vce_opts cluster(`T' `panel_id') bw(`bw')
	}
	else if "`cluster_type'"=="r" {
		* basic robust SEs
		local vce_opts robust  bw(`bw')
	}
}
***********************************

* run regressions, rolling across maturities
local idx 0
local curr_min 0
while `curr_min' <= `mat_window_final' {
	n di "`curr_min'"
	local idx = `idx' + 1
	* include maturities between curr_min and curr_max
	local curr_min = (`idx'-1)*`mat_step'
	* optionally increased window size for long maturities
	if `curr_min' < 20 {
		local curr_max = `curr_min' + `mat_window'
	}
	else {
		local curr_max = `curr_min' + `mat_long_window'
	}
	* regression shock on shock star
	`reg_cmd' `D_shock' B_shock_S0 B_shock_S1 B_shock_L0 B_shock_L1 ///
		if `mat_var'>=`curr_min' & `mat_var'<`curr_max', ///
		a(`panel_id') `vce_opts'
	* save results to vector (dropping constant)
	mat b_vec = e(b)
	mat b_vec = b_vec[1, 1..4]
	* flatten variance-covariance matrix
	mat V_mat = e(V)
	mat V_mat = V_mat[1..4, 1..4]
	mata: st_matrix("V_vec", colshape( st_matrix("V_mat"), 16) )
	* get col/row names
	local V_rownames: rownames V_mat
	local V_colnames: colnames V_mat
	local V_vec_colnames ""
	forval r=1/4 {
		forval c=1/4 {
			local Vrn: word `r' of `V_rownames'
			local Vcn: word `c' of `V_colnames'
			local V_vec_colnames `"`V_vec_colnames' "V_`Vrn'_`Vcn'" "'
		}
	}
	matrix colnames V_vec = `V_vec_colnames'
	* additional stats
	* test equality of short/long coeffs
	test B_shock_S0 = B_shock_L0
	scalar p0 = r(p)
	test B_shock_S1 = B_shock_L1
	scalar p1 = r(p)
	mat reg_stats = (`curr_min', `curr_max', e(N), e(N_clust), ///
		e(df_r), e(r2_a_within), p0, p1)
	matrix colnames reg_stats = "maturity_start" "maturity_end" ///
		"N" "N_clust" "df_r" "r2_a_within" "eq_pval0" "eq_pval1"
	* combine
	mat res = (reg_stats[1, 1...], b_vec[1, 1...], V_vec[1, 1...])
	
	* start with row of all zeros
	if `idx'==1 {
		matrix results_mat = 0 * res
	}
	* append
	matrix results_mat = (results_mat \ res)
	
}
cap drop B_shock*


* save results
if "`fname_stub'" != "" {
	preserve
	clear
	svmat results_mat, names(col)
	* save to dta file
	save ../data/rolling/phaltreg_`fname_stub'.dta, replace
	restore
}



end




